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We use Monte Carlo simulations of the 2D one component Coulomb gas on a triangular lattice, 
to study the depinning transition of a 2D vortex lattice in a commensurate periodic potential. A 
detailed finite size scaling analysis indicates this transition to be first order. No significant changes 
in behavior were found as vortex density was varied over a wide range. 
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The theory of defect mediated melting of a two dimen- 
sional (2D) solid, introduced by Kosterlitz and Thouless 
(KT) (| and developed by Nelson and Halperin (NH) @] 
and Young j|| in the late 70's, has remained a topic of 
active investigation. Numerous numerical studies have 
reported conflicting results as to whether the 2D melting 
transition is indeed a 2nd order KT type transition, or 
whether it is 1st order |Q. The most recent results have 
supported the KT scenario §. In the same paper I2| 
in which they developed this theory of 2D melting in a 
continuum, Nelson and Halperin also considered the case 
of a 2D solid in a commensurate periodic potential. In 
this case, they argued that the 2D solid would have two 
distinct phases: a "pinned" solid with long range transla- 
tional correlations at low temperatures, and a "floating" 
solid with algebraic translational correlations, similar to 
that found in the continuum, at intermediate tempera- 
tures. By mapping the problem to a vector Coulomb gas, 
they argued that this depinning transition was a 2nd or- 
der KT type transition similar to that of melting, with 
a universal discontinuous jump in the exponent rj char- 
acterizing the algebraic correlations of the floating solid 
phase. 

Despite the wide attention given to the continuum 2D 
melting problem, this depinning transition has been very 
little studied. Only recently has the floating solid phase 
been observed in numerical studies of 2D vortex lattices 
H , and related XY models Q . In this paper we present 
the first detailed finite size scaling analysis of this depin- 
ning transition. We treat the specific case of logarithmi- 
cally interacting 2D vortices. This is not only a system 
of considerable recent interest in connection with high 
temperature superconductors, but also has a unique ad- 
vantage for numerical simulations: it is an incompressible 
system. We can simulate at constant density, yet there 
will still be a sharp transition temperature, rather than 
the finite temperature interval of coexisting phases that 
one would have for other interactions (should the tran- 
sition be 1st order). We can therefore avoid the contro- 
versy, that arose in 2D melting simulations, as to whether 
one should use a constant pressure rather than a constant 
volume ensemble. For this vortex system, our results are 



quantitatively consistent with a 1st order depinning tran- 
sition. 

Our model is the one component classical Coulomb gas 
on a triangular lattice, given by the Hamiltonian, 

n = \Y i (n i -f)G ij (n j -f) . (1) 

Here rii = 0, 1 is the integer charge at site i of an L x L 
periodic triangular lattice; — / is a uniform background 
charge density; charge neutrality fixes the number of in- 
teger charges to Y^i n i = fL 2 - Gij is the 2D Coulomb 
potential for a discrete triangular lattice with periodic 
boundary conditions Q, and the sum is over all pairs 
of sites. Gij ~ In jr.; — r 3 | for distances large compared 
to the lattice spacing, but small compared to L. The 
charges n, model logarithmically interacting vortices in 
the phase of a superconducting wavefunction, / is the 
number of applied magnetic flux quanta per unit cell of 
the triangular lattice, and restricting the vortices to the 
sites of the triangular lattice models the periodic pinning 
potential Q. 

We study the Hamiltonian (|l]) using Monte Carlo (MC) 
simulations. Details of our simulation methods follow 
those described in Ref. ||. Our main results are for the 
fixed vortex density of / = 1/100. 

To decide the order of the depinning transition at T c , 
we study the finite size scaling of the energy density his- 
togram P(E) measured close to T c . For a first order tran- 
sition, one expects to find a bimodal P(E), with peaks 
that sharpen to two separated 5— functions as L — > oo; 
these give the differing energies of the two coexisting 
phases. In Fig. |l| we plot the normalized P{E) for sizes 
L = 120, 160, and 200, near T c ~ 0.00229. Analysis of 
the vortex structure function, 

5 ( k ) = 7l2E< elk ' (ri ~ r,) > ' ( 2 ) 

1 i,3 

clearly shows that the lower (upper) energy peak cor- 
responds to the pinned (floating) vortex lattice. For 
L = 120, 160, and 200 respectively, we have carried out 
0.77, 1.5 and 4.7 x 10 s MC passes though the lattice to 
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compute averages. These resulted in 2219, 1186 and 841 
hops, respectively, between the two peaks. Having a large 
number of hops is important to achieve good equilibra- 
tion of the relative weights of the two peaks (we have also 
done simulations for L = 80, 100, 140, and 180). Quali- 
tatively, we see that these peaks sharpen as L increases, 
consistent with a 1st order transition. We now make this 
observation quantitative using several different criteria. 
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FIG. 1. Histograms of energy density P(E) vs. E for 
vortex density / = 1/100, and system sizes L — 120, 160 and 
200, at temperatures near the depinning transition T c . The 
two peaks correspond to the coexisting pinned and floating 
solid phases. The peaks sharpen as L increases, consistent 
with a 1st order transition. The solid lines are fits to the 
peaks, as described in the text. 

For an infinite system, the bimodal P(E) should exist 
only precisely at T c , where the two phases coexist. For 
finite L, this coexistence region where P(E) is noticeably 
bimodal persists over a finite temperature interval AT. 
Since the relative weight of each peak is determined by 
the total free energy difference between the two phases, 
AT - L d , we expect the scaling AT ~ 1/AT = l/L d = 
1/L 2 in 2D. Using standard methods || to extrapolate 
the histograms of Fig.[l] to nearby temperatures, we de- 
fine the upper (lower) limit of the coexistence region, 



T c max (mm) > as the temperature at which the height of 
the lower (upper) peak has decreased to 1% the height of 
the upper (lower) peak. In Fig. |] we plot the results. We 
see that T c max and T c m in converge to a common value as 
L increases, and that AT ~ 1/T 2 as expected for a 1st 
order transition. 



H 



0.0026 
0.0025 
0.0024 
0.0023 
0.0022 
0.0021 
0.0020 



(a) 



60 



80 



100 120 



14000 



140 y 160 



180 200 220 




10000 



20000 L2 30000 



40000 



50000 



FIG. 2. (a) Plot of the upper and lower limits, T cma x and 
T c m i n , of the coexistence region, and the finite size depinning 
temperature T c , vs. system size L. (b) Inverse of the width 
of the coexistence region, 1/AT = l/(T cmax — T cm i n ), vs. 
L 2 . Vortex density is / = 1/100. The straight line is a least 
squares fit, and the good agreement is consistent with a 1st 
order transition. 

Next, we consider in detail how the two peaks sharpen 
as L increases. For a 1st order transition at T c , the 
width a of each peak is determined by ordinary non- 
critical finite size fluctuations giving, a ~ 1/L d l 2 = 1/L 
in 2D. The total width of the bimodal P(E), 0~t o t, ap- 
proaches a constant (this is equivalent to the familiar 



observation that the specific heat 



- Td 



L d a tot /T ~ L d 



at a 1st order transition). To verify these scaling be- 
haviors, we need to deconvolve the P(E) of Fig. |l| into 
two separate peaks, and to determine T C (L) by the cri- 
teria that the two peaks subtend equal areas ||. For 
large enough T, we expect each peak to have a Gaus- 
sian shape, and we have found this to give a good fit 
for the floating solid peak. The pinned solid peak how- 
ever lies too close to the ground state energy E ; even 
for our biggest size L = 200, the states in this peak cor- 
respond to a few discrete excitations above the ground 
state. We therefore fit the pinned solid peak to the ad 
hoc form, [E - T )exp[-(T - Ti) 2 /2er 2 ], where E x and 
a are fitting parameters. Using data from only the far 
side of each peak, we fit P(E) to the sum of a Gaus- 
sian and the modified Gaussian above, to get the fitted 
curves shown in Fig. [l| From these fits we determine the 
probabilities that a state with a given value of E belongs 
to the pinned phase, the floating phase, or neither (the 
later being the "transition" states, consisting of large do- 
mains of one phase in a background of the other, which 
give rise to the transitions between the two phases). We 
then go though the ensemble of states that enter into our 
averages, and probabilistically assign each state to the 
pinned phase, the floating phase, or neither. With this 
deconvolution, we then extrapolate to the T C (L) which 
gives equal weight to the two phases, and we can then 
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compute the widths a^ oa ±, 0-pin, and er t ot at this tem- 
perature. Our result for T C {L) is shown in Fig.|^a. Our 
results for the widths are shown in Fig.||. We find er t ot ~ 
constant, and <7fl oa t, Cpin ~ 1/L, as expected for a first 
order transition pi] . We can also use this decomposition 
to compute the entropy jump per vortex at the transi- 
tion, As = ((-E)float - {E)p in )/(fT c ). Our results, shown 
in Fig.|^, give a As(L) that saturates to a constant as L 
increases, again as expected for a 1st order transition. 



V(T+) = 4/ , 



(4) 
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FIG. 3. Inverse of the total width a to t of the energy his- 
togram P(E), and inverse of the widths, a P i n and erfloat, of 
the pinned and floating solid peaks, vs. system size L. Vortex 
density is / = 1/100. The straight lines are least squares fits. 
The observed scaling l/a to t ~ constant, l/opin, 1/cfloat ~ L 
is consistent with a 1st order transition. 



where / is the vortex density. From Eq.(^|) we see that 
5(Ki)/i 2 vs. L plotted on a log-log scale should yield 
a straight line of negative slope rj(T). In Fig. Ha we 
show such plots for several values of T near T c , using 
S(Ki) averaged over all configurations encountered in 
the simulation. The straight lines are a least squares fit 
to the assumed algebraic form. We see that the fit is not 
very good, except perhaps at the highest temperatures. 
If, however, we average iS(Ki) separately over only the 
states in the floating solid phase, and over only the states 
in the pinned solid phase, we see the expected behavior 
as shown in Fig.J^b. In the pinned phase, the curves 
saturate to a finite value as L increases, reflecting the 
long range order of this phase. In the floating phase, the 
curves give good straight line fits, showing algebraically 
decaying correlations. From these fits, we extract the 
exponent 77 of the floating phase, which we plot vs. T 
in Fig.^j. We see that r/(T+) is close to, but noticeably 
above, the NH prediction of Eq.([|). Since the NH map- 
ping onto a vector Coulomb gas inverts the temperature 
scale, this is consistent with a 1st order transition pre- 
empting the NH defect unbinding transition, with Eq. (|4|) 
serving as a lower bound on rj(T). 
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FIG. 4. Entropy jump per vortex, 

As — ({i?)fl at — (-E)pi n )/(/T c ), at the depinning transition 
vs. system size L. The vortex density is / = 1/100. The 
dashed line is a guide to the eye only. 

Next we consider the translational correlations. In 
the floating solid phase correlations are algebraic 0, 
(exp[iK • (rj — r^)] ~ |r$ — Yj \~ T > K ( T \ with K a recip- 
rocal lattice vector, and t\k = ??|K| 2 /|Ki| 2 , with Ki the 
smallest non-zero K. Substituting this into the structure 
function, Eq.(|2|) results in, 



S(K) 
L 2 



-Vk 



-(7,lnL)|K| 2 /|Ki| 2 



(3) 



According to the NH theory jj, r\ should take a discon- 
tinuous jump from zero in the pinned solid to a universal 
finite value at T c in the floating solid. For an incom- 
pressible system such as ours, in which the bulk modulus 
A = 00, the predicted jump is, 
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FIG. 5. Log-log plot of S'(Ki) vs. L, for several values 
of temperature near T c . (a) Shows S'(Ki) averaged over all 
configurations; (b) shows S'(Ki) averaged separately over only 
the pinned and only the floating configurations. From top to 
bottom the temperatures are, T = 0.00222, 0.00226, 0.00230, 
0.00234 and 0.00236. 
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FIG. 6. Translational correlation exponent rj(T) of the 
floating lattice phase, as obtained from the fits in Fig. |E| The 
vortex density is / = 1/100. rj(T c ) remains above the theo- 
retical lower limit 4/ of the Nelson Halperin theory. 

The results above are all consistent with a 1st order 
transition rather than the 2nd order NH prediction. How- 
ever one can question whether this is an artifact of the 
particular density / = 1/100 that we have used. Hattel 
and Wheatley jll[] in particular have argued that, as / 
decreases in the vortex system, the core energy of the de- 
fects that lead to depinning increases, and hence even if 
depinning is 1st order at some /, it must ultimately be- 
come 2nd order as / decreases. To check this prediction 
we have carried out simulations for a variety of densi- 
ties / = 49, 64, 100, 196, 400. Rather than do finite 
size scaling for each case, we study only systems with a 
fixed total number of vortices N v , corresponding to sys- 
tem sizes L = y/N v /f. We choose N v = 144 because our 
finite size studies of / = 1/100 indicated this to be suf- 
ficiently large to be within the asymptotic large L limit. 
Analyzing the data from these simulations in exactly the 
same manner as described above, we show in Fig. ^ our 
result for the entropy jump per vortex at depinning, As, 
as a function of density /. We see that As is only weakly 
dependent on /, extrapolating to a finite As ~ 0.60 as 
/-0. 
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FIG. 7. Entropy jump per vortex, 

As = ((E) fl oa t — (E) p i n )/(fT c ), at the depinning transition 
vs. vortex density /, for systems with a fixed number of vor- 
tices N v = fL 2 = 144. The dashed line is a least squares fit 
showing only a weak decrease in As as / decreases. 



is slightly larger than the more correct determination in 
Fig. H However the main point of Fig.|| is to observe that 
the estimates rj' all collapse to a common curve as the 
density / decreases. There is no evidence for a continuing 
decrease in the value of r)(T+) as / decreases, as might 
be the case if the 1st order transition was weakening and 
the NH prediction of Eq. (H) was approached more closely. 
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FIG. 8. Estimated exponent r)' tot /f, r)' pin /f, and 7?fl oat //, 
as determined from S(K) averaged over all states, the pinned 
solid states, and the floating solid states, respectively, vs. 
T/T c (f) for different vortex densities /. The data collapse 
to a common curve as / decreases. 

To conclude, a finite size scaling analysis for the den- 
sity / = 1/100 is strongly consistent with the depinning 
transition being 1st order. We further find no evidence 
that this 1st order transition changes in any significant 
way if the density / is decreased. We cannot rule out 
the possibility of different behavior for softer interactions 
than the logarithm considered here. 

We would like to thank D. R. Nelson and M. Franz for 
very helpful conversations. This work was supported by 
the Engineering Research Program of the Office of Basic 
Energy Sciences at the Department of Energy, grant DE- 
FG02-89ER14017. 



Finally, we consider the translational correlations for 
the different /. Since we have simulated only for fixed 
values of fL 2 , we cannot extract r\ by scaling with L as 
in Fig. ^|b. However we can estimate rj as follows. From 
Eq.(||) we expect a Gaussian envelope for the peaks S(K) 
as a function of |K|; the width gives a measure of r\. This 
envelope is not purely Gaussian; a correction exists due 
to a K dependence of a prefactor in the L~ VK scaling of 
Eq.(^|). Making a simple approximate correction for this 
effect, as done in Ref. [||, we plot our estimate for the 
correlation exponent, rf, in Fig.|^. We show results using 
S(K) averaged over all states, only the floating states, 
and only the pinned states, plotting the data in scaled 
units of rj'/f vs. T/T c . For / = 1/100, our estimate ?7g oat 
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